function [F peakLocation P countPeaks] = peak31finder(X,T,r,sigma,n)
P=zeros(n,3);
F=zeros(r+1,r+1);
peakLocation=zeros(r+1,r+1);
%%derivativeLocation=zeros(r+1,r+1);
h=pi/(2*r);
countPeaks=0;
for i=0:r
    theta=i*h;
    for j=0:r
        phi=j*h;
        F(i+1,j+1)=f31(theta,phi,sigma,T,X);
    end
end

for j=2:r
    if((F(1,j)>F(2,j-1))&&(F(1,j)>F(2,j))&&(F(1,j)>F(2,j+1))&&(F(1,j)>F(1,j-1))&&(F(1,j)>F(1,j+1)))
        peakLocation(1,j)=1;
        countPeaks=countPeaks+1;
        P(countPeaks,:)=[0 0 1];
    end
    if((F(r+1,j)>F(r,j-1))&&(F(r+1,j)>F(r,j))&&(F(r+1,j)>F(r,j+1))&&(F(r+1,j)>F(r+1,j-1))&&(F(r+1,j)>F(r+1,j+1)))
        peakLocation(r+1,j)=1;
        countPeaks=countPeaks+1;
        P(countPeaks,:)=[cos((j-1)*h) sin((j-1)*h) 0];
    end
end

for i=2:r
    if((F(i,1)>F(i,2))&&(F(i,1)>F(i-1,1))&&(F(i,1)>F(i+1,1))&&(F(i,1)>F(i-1,2))&&(F(i,1)>F(i+1,2)))
        peakLocation(i,1)=1;
        countPeaks=countPeaks+1;
        P(countPeaks,:)=[sin((i-1)*h) 0 cos((i-1)*h)];
    end
    if((F(i,r+1)>F(i,r))&&(F(i,r+1)>F(i-1,r+1))&&(F(i,r+1)>F(i+1,r+1))&&(F(i,r+1)>F(i-1,r))&&(F(i,r+1)>F(i+1,r)))
        peakLocation(i,r+1)=1;
        countPeaks=countPeaks+1;
        P(countPeaks,:)=[0 sin((i-1)*h) cos((i-1)*h)];
    end
end
if((F(1,1)>F(1,2))&&(F(1,1)>F(2,1))&&(F(1,1)>F(2,2)))
    peakLocation(1,1)=1;
    countPeaks=countPeaks+1;
    P(countPeaks,:)=[0 0 1];
end

if((F(1,r+1)>F(1,r))&&(F(1,r+1)>F(2,r))&&(F(1,r+1)>F(2,r+1)))
    peakLocation(1,r+1)=1;
    countPeaks=countPeaks+1;
    P(countPeaks,:)=[0 0 1];
end

if((F(r+1,1)>F(r,1))&&(F(r+1,1)>F(r+1,2))&&(F(r+1,1)>F(r,2)))
    peakLocation(r+1,1)=1;
    countPeaks=countPeaks+1;
    P(countPeaks,:)=[1 0 0];
end

if((F(r+1,r+1)>F(r,r))&&(F(r+1,r+1)>F(r,r+1))&&(F(r+1,r+1)>F(r+1,r)))
    peakLocation(r+1,r+1)=1;
    countPeaks=countPeaks+1;
     P(countPeaks,:)=[0 1 0];
end

for i=2:r
    for j=2:r
        if((F(i,j)>F(i-1,j-1))&&(F(i,j)>F(i-1,j))&&(F(i,j)>F(i-1,j+1))&&(F(i,j)>F(i+1,j-1))&&(F(i,j)>F(i+1,j))&&(F(i,j)>F(i+1,j+1))&&(F(i,j)>F(i,j-1))&&(F(i,j)>F(i,j+1)))
            peakLocation(i,j)=1;
            countPeaks=countPeaks+1;
            P(countPeaks,:)=[cos((j-1)*h)*sin((i-1)*h) sin((j-1)*h)*sin((i-1)*h) cos((i-1)*h)];
        end
    end
end